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Abstract. - We perform Transition matrix Monte Carlo simulations to evaluate the entropy 
of rhombus tilings with fixed polygonal boundaries and 2D-fold rotational symmetry. We 
estimate the large-size limit of this entropy for D — 4 to 10. We confirm analytic predictions 
of N. Destainville et ai, J. Stat. Phys. 120, 799 (2005) and M. Widom et al, J. Stat. Phys. 
120, 837 (2005), in particular that the large size and large D limits commute, and that entropy 
becomes insensible to size, phason strain and boundary conditions at large D. We are able to 
infer finite D and finite size scalings of entropy. We also show that phason elastic constants 
can be estimated for any D by measuring the relevant perpendicular space fluctuations. 



Random tiling models [1,2] have been intensively studied since the discovery of quasicrys- 
tals in 1984, because they are good paradigmatic models of quasicrystals. These metallic 
compounds exhibit exotic symmetries {e.g. icosahedral) which are classically forbidden by 
crystallographic rules. This is accounted for by the existence of underlying quasiperiodic or 
random tilings. When tiles are decorated in some manner by atoms, these tilings become ex- 
cellent candidates for modeling real quasycristalline compounds [3] . As compared to perfectly 
quasiperiodic tilings, some specific degrees of freedom, the so-called phason flips, are activated 
in random tilings, giving access to a large number of microscopic configurations. The number 
of configurations grows exponentially with the system size in contrast to perfectly quasiperi- 
odic tilings where it only grows polynomially. The resulting configurational entropy lowers 
the free energy as compared to competing crystalline phases. Despite their random character, 
random tilings still display the required macroscopic point symmetries in their Fourier spectra. 
They are as good candidates as perfectly quasiperiodic tilings for modeling quasicrystals [1,2]. 
The statistical mechanics of random tilings is of central interest for quasi-crystallography. But 
the calculation of the thermodynamical observables of interest in random tiling models has 
turned out to be a formidable task. Even the calculation of configurational entropy in the case 
of maximally random tilings where all tilings have the same energy is a notoriously difficult 
problem. Very few models are exactly solvable [4-9], and a large majority of calculations rely 
on numerical simulations. In Refs. [10,11], an original analytic mean-field theory for plane 
rhombus tilings with 2D-fold symmetry, in the large D limit, was proposed. Its ultimate 
goal is to derive valuable results on finite D tilings by estimating finite D corrections to the 
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infinite D limit. The present study intends to support numerically the analytic findings of this 
promising approach. We propose a calculation of the entropy per tile of rhombus tilings based 
on a Transition matrix Monte Carlo (TMMC) technique [12]. We get estimates of the large 
size limit of the entropy for D < 10 (TableHJ. Furthermore, we calculate the same entropy for 
large D tilings filling 2_D-gons of finite side length p. We confirm the theoretical predictions 
of Refs. [10, 11]: large D and large p limits commute. We also explore the large D behavior 
of a phason elastic constant by measuring the relevant perpendicular space fluctuations. We 
show that it decays like 1/D. 

We study tilings that fill centro-symmetric 2£>-gons (Fig.^). Such fixed boundaries have 
a strong influence on tilings at finite D (see [11, 12] and references therein). Their most 
spectacular consequence is that the entropy per tile of fixed-boundary tilings is strictly smaller 
than its free- or periodic-boundary counterpart. It was also one of the goals of Ref. [11], where 
these questions are thoroughly discussed, to explore what this difference between fixed and 
free boundary conditions becomes in the large D limit. It was demonstrated there that, in this 
limit, all boundary conditions become equivalent. More generally, it was proved that under 
very weak conditions, the entropy per tile of large D tilings is independent of their size, shape, 
tile fractions and boundary conditions. The "universal" entropy per tile was estimated by a 
mean-field approach and numerically calculated : ~ 0.5676 ± 0.0001. It was also argued 
that some phason elastic constants associated with given strain modes vanish at large D. Our 
goal here is to explore numerically these predictions in a larger variety of situations. 

The tiles considered in this paper are lozenges of unitary side length and of angles multiple 
of ir/D. Their edges are collinear to the D vectors e|, a = 0, . . . , D — 1 which make angles 
cm I D with an arbitrary direction. In the celebrated cut-and-project framework (see [1]), rhom- 
bic tiles are considered as the projections onto the 2-dimensional space of the 2-dimensional 
faces of a hypercube in a space of dimension D denoted by H. Conversely, any tiling in the 
bidimcnsional space can be "lifted" to a bidimensional continuous membrane embedded in 
the Z D lattice. When this membrane is projected back into the 2-dimensional space, the 
facets of the Z?-dimensional lattice project precisely onto the rhombic tiles. The normalized 
basis vectors of H are denoted by e a . The e& are the normalized projections of the e Q [10]. 
The physical 2-dimensional space is usually called the parallel space and it is denoted here 
by The (D — 2)-dimensional space in H perpendicular to is called the perpendicular 
space. It is denoted by E 1 - . The normalized projections of the e a in E 1 - are the e^. The 
2D-gons filled by such tiles in the present publication have edges parallel to the e|. Their 
integral side lengths are denoted by p a . When all p a are equal to the same p, the tilings 
are said here to be diagonal, of side p. Local moves called single flips are usually used to 
explore configuration spaces of rhombus tilings. In dimension 2, they consist of the exchange 
of three tiles that locally fill a small elementary hexagon (see Fig.^ left). Such flips connect 
configuration spaces in dimension 2. 

Transition matrix Monte Carlo technique. - The previous estimates of the large D en- 
tropy [11] relied on an efficient numerical approach in the p = 1 case. It cannot, however, 
be applied to values of p > 1 of interest in the present paper, without inextricable technical 
complications. We chose to apply here a Transition matrix Monte Carlo technique [12]. It is a 
variant of the transition matrix method [13]. It uses a standard Metropolis Monte Carlo sam- 
pling to construct a numerical approximation to the transition matrix. The density of states 
is calculated from this transition matrix and provides the total number of configurations. The 
energy E used here in order to implement the Monte Carlo algorithm has no particular physical 
meaning as far as quasicrystals are concerned. It is related to the structure of the configura- 
tion space in terms of partial ordered sets (posets) theory [14]. It involves the description of 
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rhombus tilings with polygonal boundaries by means of generalized partitions, itself related 
to the cut-and-project technique. Its complete construction is detailed in Refs. [12] and [14] 
(where it is equal to the so-called "rank function"). There is a unique tiling t m i n (resp. t max ) 
with minimal energy E m i n (resp. maximal energy E max ) (see Fig. ^ bottom). The energy 
measures the distance from a given tiling t to t m in- A single flip as described above increases 
or decreases the energy by a single unit, according to whether it increases or decreases the dis- 
tance from t m i n . The structure of graded poset [14] ensures that the energy variation between 
any two tilings does not depend on the choice of the path of flips between these tilings. We 
find bellow that the energy of most tilings is close to E mi( i = (E min + E max )/2. The density 
of states is symmetric with respect to E mi d, because of the geometric centro-symmetry of the 
problem. 




Fig. 1 - Several examples of tilings of side p — 5 with D = 4. (a,b): two typical random tilings 
differing by a single flip in the emphasized hexagon. (c,d): Tilings with minimal energy E min (c) and 
maximal energy E max (d). Tilings (a) and (b) have energies close to E mid = (E min + E max )/2. 

In the Metropolis Monte Carlo algorithm, at temperature T > 0, a randomly chosen vertex 
is tried to be flipped if it is flippable. This change is accepted if it lowers the energy. If by 
contrast the change raises the energy, it is accepted with probability exp(— AE/T), where AE 
is the energy variation if the change was accepted. At low T, tilings are close to t m i n . At 
large T, they are in the maximum entropy region, close to E m id- It will also be useful to take 
negative values of T after adapting suitably the previous transition probabilities. Now we use 
the fact that the energy variation associated with a single flip cannot be larger than ±1. As 
in Ref. [12], we denote by n±(t) the number of tilings that can be reached from a given tiling 
t by single upwards or downwards flips. For any given tiling t, the numerical calculation of 
n±(t) is easy and exact. If Ny is the number of tiling vertices, then n+(t) + n_(t) < Ny. 
Remind that Ny is tiling-independent. In particular, it is unchanged through a single flip. 

Now we can define the transition matrix u>(E,E') of dimension E max — E rn i n + 1. It is 
equal to everywhere except on its diagonal and two off-diagonals: 



where P(E) is the set of tilings of energy E and W(E) is their number. In addition u>(E, E) = 
1 — U)(E,E— 1) — u>(E, E+l). The matrix u> is not known and it is the goal of our Monte Carlo 
algorithm to estimate it by sampling the energy levels E and estimating the exact mean JIJ 
by an approximate average on sampled configurations. The density of states can then be 
extracted as follows: lu + (E)W(E) = oj-{E + 1)W{E + 1) because the total number of forward 
flips from energy E to energy E + 1 is exactly equal to the total number of backward flips from 
energy E + l to energy E. This equation reads W(E + 1) = uj + (E)W(E)/lu-(E + 1). It allows 
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us to iteratively extract the W(E) using uniqueness of the ground state t m i n , W(E m i n ) = 1. 
Finally, the total number of tilings is Z = J^e W(E), the configurational entropy is S = InZ, 
and the entropy per tile is 075 (p) = S/Nt, Nt the number of tiles. To ensure that all energy 
levels are (almost) uniformly visited, we perform sweeps over temperature, for both positive 
and negative temperatures. The minimum and maximum temperatures T m i n and T max must 
be suitably adjusted to ensure that this distribution is as uniform as possible in order to sample 
correctly all energy levels. We find that the density of states W(E) is nearly Gaussian, at 
least near its maximum at E m id [12]. 

Numerical results. We first focus on the diagonal case where p a = p. Tilings are 
initialized at minimum energy £ ml „ which is chosen to be equal to 0. The maximum energy 
can be calculated a priori [12,14]: E max = ( 3 ) p 3 . We accumulate statistics n±(t) to estimate 
the transition matrix. The degree to which symmetry of S(E) with respect to E m id is broken 
serves as an indicator of errors accumulated during the iterative calculation. The residual 
R = In W(E max ) ) that should ideally be equal to 0, reflects the cumulative errors in W(E). It 
can be demonstrated [12] that R is a good indicator of the uncertainty on S at the end of the 
calculation. It will provide error bars on o~d{p) below. In all cases, error bars will be smaller 
than 10~ 4 (except for D > 50 where they will be smaller than 5.1CP 4 ). 

We begin with the infinite p limit of finite D entropies of diagonal tilings with fixed 
polygonal boundaries: Tfrj = linip^oo Wd(p)- For D — 4— 10, we extrapolate the large p limit 
by fitting the numerical data obtained via the previous algorithm. We choose a second order 
fit in I /p. Contrary to the three-dimensional case [12], we did not observe that logarithmic 
corrections improved the quality of the fit. Fig. (left) displays <td(p) m function of 1/p and 
the second order fits. The extracted large p limits are listed in Tab. |U We applied the same 
procedure to the exactly solvable D = 3 case (see [1]) and we found that the method provides 
the exact expected value with a relative precision smaller than 1%. Therefore we anticipate 
that the relative errors for larger values of D are also smaller than 1%. 
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Fig. 2 - Left: Entropies per tile Wd{p) vs 1/p for D — 4,5,6,7,8,9 and 10 from bottom to top 
(symbols) and second order fits (lines). Error bars are much smaller than symbol sizes because all 
errors are smaller than 5 10 -4 . Right: Entropies per tile vs 1/D in various situations (symbols; 
see text) and their fits at order 2 in 1/D (lines). All plots converge towards the "universal" value 
(Too ~ 0.5676 within our uncertainties. 



To test the predictions stated above, we also extrapolate the limit aoo — lim£>^oo 5\d = 
0.563±0.006 by fitting the finite D values at order 2 in 1/D [11]. The error bar also corresponds 
to a relative error of 1%. We get Tjoo = within error bars. This is a first confirmation of 
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Table I - Large size limits of the entropy per tile. Relative errors are smaller than 1% (see text). 



D 


4 


5 


6 


7 


8 


9 


10 


a(D) 


0.360 


0.410 


0.441 


0.461 


0.476 


0.487 


0.496 



the predictions under interest. To verify them further, we also calculate the large D limit for 
finite p. The case p = 1 was solely considered in Ref. [11]. We were able to perform more 
intensive numerical calculations here. The entropy is extrapolated from finite D entropies for 
D = 20 to 70 for p = 2, and D = 20 to 50 for p = 3. We find consistently a^ip = 2) = 
&oo(p — 3) = 0.568 ± 0.006. Then we consider boundary strains [11], with high frequency, 
(jPa) = (l,f>, .. . , l,p), with p = 2, 3, as well as low frequency, (p a ) = (1, . . . , l,p, . . . ,p), with 
D/2 sides of length 1 (resp. of length p = 2,3). Extrapolating data for 20 < D < 50, we 
also find <?oo(l, 2, 1,2,...) = CToo(1, 3, 1,3,...) = 0.569 ± 0.006 at high frequency strain and 
(7oo(l, . . . , 2, . . .) = CToo(l, . . . , 3 . . .) = 0.568 ± 0.006 at low frequency. These calculations, 
summarized in Fig. (right), also support the predictions of the "universality" of the large D 
limit [10, 11]. They also confirm the commutation of the large p and large D limits. 

Now we examine finite p and finite D scalings. A priori, if there are not logarithmic 
corrections, the Taylor expansion of a(D,p) at order 2 in powers of 1/p and 1/D reads 

,„ . a b c d e 

a(D,p)~a 00 + - + - + — + — + —. 2 
p D p z p D D z 

Since a(oo,p) = a Xl a = c = 0. Moreover, at fixed D = Do and at first order in p, 
7t(Do,p) ~ (cfoo + tj^ + yj?) + 7^7 \- From the fit of the data displayed in Tab. |U we get 
b ~ —0.57 and e«-l. It is not possible to get a reliable estimate of d by fitting the available 
small Dq data. Indeed, it appears that the slopes dj Dq do not have reached their asymptotic 
regime at D = 10 because they are still increasing with Dq. However, fitting these available 
slopes at the second order in 1/Dq provides the rough estimate d w 0.2. In addition, at fixed 
p = Po and at first order in 1/D, Eq. J3J becomes a(D,p ) ~ CTco + (6 + d/p )/D. Fitting 
the data for po = 2 and po = 3 gives respectively b + d/2 ~ —0.31 and b + d/3 ~ —0.41, 
from which we get b ~ —0.61 and d ~ 0.60. These values are compatible with the previous 
ones. In particular, both values of b match very well, whereas they are extracted from the 
numerical data in a different way: in the first case, we take the infinite p limit first, and in the 
second case, the infinite D limit first. This confirms a posteriori the general expression J2J. 
In particular the two limits commute; There are probably no logarithmic corrections in (|2"|). 

Phason elastic constants. - Now we use the previous algorithm to estimate phason elastic 
constants. Even though we are working on fixed-boundary tilings, our goal is to measure elastic 
constants associated with /ree-boundary ones. We use the notations and the conventions of 
Ref. [10]. In particular, the normalizations must be precisely prescribed: we define 3 scalar 
products in H, £" and denoted by (.|.)y where V is the space under consideration. If 
x,y e £-S (x\y) £ ± = s\(x\y) H with s\ = D/{D-2) [10]. In a similar way, sjj = D/2. With 

these conventions, we get ||e Q ||ij = ||e^|| £ i = ||ei|| £ n = 1, as desired. 

We also employ here a useful variational principle (see [12] and references therein) to write 
the entropy in a practical field-theoretic manner. After coarse-graining [1], large size tilings 
are represented by regular height functions <j> : R 2 — > R D ~ 2 . Such a "macroscopic" variable 
represents a large number N(<p) of tilings, whose membrane representation in Z D is equal, 
after coarse-graining, to the graph of <p. Their micro-canonical entropy \og(N(<f))) is equal to 
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a functional S[tf>], which is the integral over the tiled domain S of a function ct(V</>) of the 
gradients of (f> [1]. The maximization of S[cj>] gives a (unique) function <fi^ which represents 
the dominant tilings in the statistical ensemble under consideration. The random tiling model 
states that cr(V</>) has a unique maximum, corresponding to tile fractions maximizing (in the 
present case) the 2£>-fold rotational symmetry. The orientation of the real space £H is chosen 
so that the gradients are zero at this maximum, and the model states that the entropy density 
has a quadratic behavior near this maximum: cr(V0) = o- m ax — \ • K • V</> + o(|V</>| 2 ). 
By analogy with an usual elastic theory, the tensor K is called the tensor of phason elastic 
constants. In the basis of the perpendicular space E 1 - where K is diagonal, the previous 
quadratic form becomes cr(V0) = (T max — | Y^i=o K i(^0i) 2 + °(\^4>\ 2 )- We denote by Oi the 
normalized (in £-"-) eigenvector of K corresponding to the eigenvalue «j. We set 

Oi[<j>] = J<f>(r) ■ Oi al 2 r. (3) 

Oi [<p] is the mean height of <fi in the direction Oi . As written above, the total entropy of a 
function <p is S[tp] = J s cr(V0(r)) d 2 r. To calculate the entropy S(Oi) of the tilings of £ whose 
mean height is Oi, we introduce the Lagrange multiplier \ and we now maximize Gi[<p] = 
S{4>] + \iOi{4>}. Without loss of generality we focus on Go- Then differentiating Go with 
respect to 0, we get, at first order in cj> and its gradients, the equations |^ = n A(j)o + A = 
and ^ a = kqA^ = for i > 0, together with the condition that <fi matches the boundary 
conditions induced by the tiling ones. Let 4>^> be the function maximizing S[4>], which also 
satisfies Oo[</>^] = 0, and <)P^ = <j> — (f)^ . By construction, 0^ = at the boundary of 
E. The only non-zero coordinate of cj)^ is 4>^\ It satisfies A^q 1 ^ = — Xq/k . At the large D 
limit, S tends towards a disk of radius R — pD/ir [10, 11]. Taking into account the condition 

4> { o ] = on this circle, we get ^ 1} (r) = -j^(r 2 - R 2 ). Thus O o [0 (1) ] = -Jr^ 4 and 



Ao = — ^ K °^° ■ One finally gets the relation (generalized to any i): 
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S(Oi) = + (D] = 5(0) - \s^Ki (4) 

for a diagonal tiling of side p with (large) 2D- fold symmetry. Therefore at T = oo, the 
equilibrium fluctuations of Oi satisfy {Of) — ^ n s K . because of the Gaussian character of the 
number exp(S(0,)) of tilings with mean height Oi. The measure of (Of) gives access to m. 

Before giving our numerical results, we discuss the validity of this calculation. It relies 
on the quadratic approximation above, which is itself valid if \7<p^ ■ K • X7<p^ <c a max = 
O(l). Now \\7cj)^\ 2 = 0(1/ D) in finite D tilings and phason elastic constants are bounded 
[10, 11]. Thus the quadratic approximation is valid at large D and we get asymptotically exact 
estimates of the Hi. But at finite D, this approach cannot be expected to provide accurate 
results. For example, if D = 3 one calculates that IV^ ^ 2 can be as large as 4 in large 
macroscopic regions near the boundary, and the previous argument collapses. We exemplify 
this procedure on Z = sJD D ~ 2 Y^a=a (~l) ae o= [!]■ O ne checks that H^H^x = 1. It is an 
eigenvector of K if and only if D is even because it belongs to a one-dimensional irreducible 
representation (irrep) in H of the symmetry group of the 2_D-gon, CiDv This irrep exists only 
if D is even. Z is a discrete component of E 1 - since for any a, {0 z \e-^} £ ± = ±.\j\jD — 2 and 
any tiling vertex is represented in E 1 - by a point of coordinate multiple of \j\jD — 2 along Z . 
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Table II - Large p estimates of the phason elastic constant k z ■ 
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11 




0.76 


0.59 


0.45 


0.38 



Tab. Illl gives the so-obtained values of k z in the large p limit for D = 5,7, 9, 11. The limits 
are estimated by a fit of k z at order 2 in 1/p. They clearly indicate that k z oc 1/D for the 
values considered here. Thus our conclusion is that n x vanishes at large D, proportionally to 
1/D, as predicted in Ref. [10]. However the elastic constants studied there were different from 
k z since they were associated with fluctuations of the de Bruijn lines separations, and existed 
whatever the parity of D. This result is in agreement with the universality of the large D 
entropy per tile: since |V</>'°'| 2 = 0(1/D), the corrections to <j x vanish at large D like 1/D 2 . 

We have mentioned that this method cannot be expected to give reliable results at finite 
D and it is indeed what we observed for D = 3. However, for D = 5, we already find, with 
the conventions of normalization of [1], k z ~ 0.26. This value is close to the expected one 
k z — 0.29 [1]. In principle it is possible to apply the same procedure to other phason elastic 
constants But determining the corresponding 9i is a tedious task, requiring a complete 
description of E 1 - in terms of irreps of C2DV This goes beyond the scope of the present paper. 

Conclusion. - The computational power provided by the Transition matrix Monte Carlo 
technique allowed us to investigate random tilings with 2Z?-fold rotational symmetry in their 
large size, large D limit. The numerical entropies obtained strongly support the analytic 
predictions of previous publications [10, 11] and go into the direction of a large D "universal" 
entropy independent of size, shape, tile fractions and boundary conditions. This work shows 
that the knowledge of a few adjustable parameters that can be estimated by finite D and 
finite p fits is sufficient to get a good estimate of finite D entropies of physical interest (see 
Eq. (J5J). For example, the fact that the numerical value of |e| is of the order of 1 shows 
that the expansion of the finite D, infinite p, fixed-boundary entropies at the first order in 
1/D is accurate within few percents as soon as D > 4. By contrast, the direct transposition 
of the present results to the more physically relevant, finite D free-boundary tilings, is less 
immediate since it requires the numerical knowledge of phason elastic constants. 

We have demonstrated that the estimation of these phason elastic constants for free- 
boundary, finite D tilings is feasible with a good accuracy. However, we did not tackle this 
task in the general case but we have demonstrated its feasibility in a general framework. In 
the case studied here, the phason elastic constant k z falls off like 1/D, in agreement with 
earlier predictions [10]. 

This work also illustrates that even though polygonal boundaries are not physical [12], 
fixed-boundary tilings are adapted to get valuable information on their free- or periodic- 
boundary counterpart. They present the great advantage of being conveniently coded and 
manipulated in a computer memory, in particular in relation with the partition techniques 
used in several publications (see [11,12,14] and references therein). 

* * * 
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salute the participation of Suriya Kothandaraman to the present work as part as his Master's 
project. 
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